Compact phases of polymers with hydrogen bonding 
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We propose an off-lattice model for a self-avoiding homopolymer chain with two different com- 
peting attractive interactions, mimicking the hydrophobic effect and the hydrogen bond formation 
respectively. By means of Monte Carlo simulations, we are able to trace out the complete phase 
diagram for different values of the relative strength of the two competing interactions. For strong 
enough hydrogen bonding, the ground state is a helical conformation, whereas with decreasing hy- 
drogen bonding strength, helices get eventually destabilized at low temperature in favor of more 
compact conformations resembling /3-sheets appearing in native structures of proteins. For weaker 
hydrogen bonding helices are not thermodynamically relevant anymore. 
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The collapse of a self-avoiding flexible polymer chain in 
a "bad" solvent has been studied for many years |l| . Fol- 
lowing de Gennes seminal work on showing the intimate 
connection of polymer collapse with tricritical systems 
, most of the theoretical effort has been concerned with 
the universal features of the #-point, the thermodynamic 
second order transition between the swollen and the com- 
pact phase |3|, this last phase beingusually regarded as 
a structureless liquid globule phase [l| . The possibility of 
more complex behavior in the compact phase has been in- 
vestigated only recently, revealing the existence, at lower 
temperatures than the collapse gas-to-liquid transition, 
of a liquid-to-solid and a solid-to-solid transition Q . 

On the other hand, protein molecules undergo simi- 
lar transitions between denatured, molten globule, and 
native states, which are solid-like structures with a well 
defined three-dimensional conformation 0. The main 
driving force of protein collapse is believed to be the hy- 
drophobic effect, which shields most of the non-polar side 
chains in the core of the native protein structure from wa- 
ter . This could indeed be grossly described as a "bad" 
solvent effect. Yet, native structures of proteins are very 
peculiar, when compared to typical compact conforma- 
tions of self-avoiding polymer chains. The benchmark of 
protein nativeness is perhaps the ubiquitous presence of 
highly ordered local motifs, called secondary structures, 
known to be stabilized by hydrogen bonding Q ■ 

In this Letter, we propose a minimal off-lattice ho- 
mopolymer model, where a usual two-body isotropic at- 
tractive interaction -mimicking the hydrophobic effect- 
is competing with a directed attractive interaction mim- 
icking the angular dependence of hydrogen bonding ||. 
We consider a chain of N beads at positions in the 
three-dimensional continuum space R 3 , with 1 < i < N. 
The chain constraint is enforced strictly, by keeping the 



distance between consecutive beads along the chain con- 
stant and unitary, — = 1, for 2 < i < N, while 
no other constraint is considered. We model the hy- 
drophobic effect @ by considering a pair-wise attrac- 
tive square well potential with a hard wall, Eh p = 
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where a is the hard-core radius of each bead, and A is 
the range of the attractive interaction. In the following 
we will always consider the case er = 1, A = 1.5, as in Q. 

In order to model hydrogen bonding, we need to 
break isotropy and favor a preferred direction between 
the two 'hydrogen-bonded' beads. We use the same 
type of directed interaction proposed by Chen and 
Kemp [JJj, so that the two planes, each containing 
one of the two hydrogen-bonded beads and its near- 
est neighbors along the chain, will both be prefe ribly 
orthogonal to the contact vector between them |ll| : 
E h b = T,2<i<j+i<N V hb (n - fj,Ui,Uj), where Uj = 
(fi+i - n) x {fi - n-i), and: 



V hb (f, u t ,u 3 ) = 0.5 (|f • fiiP + I* ' fyl" 1 ) V h p (H) , (2) 

where • denotes normalized vectors. The directionality 
degree of hydrogen bonding is controlled by the exponent 
to; a large value corresponds to a strong 'directionality'. 
We have mainly studied the to = 12 case, since lower 
values of to do not favor protein-like secondary structures 
in this parametrization. 

The interplay between hydrophobic collapse and hy- 
drogen bonding is controlled by the relative strength a 
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between the two interactions when the following total 
Hamiltonian is considered: 



~H a — Ehp + aEhb ■ (3) 

Whenever two beads come into contact, i.e. their mu- 
tual distance falls within the well, they always gain a neg- 
ative unitary energy contribution from E} lp . A further 
negative contribution could come from E^b, depending 
on how well the hydrogen bond is formed between them, 
ranging from 0, in the worst case, to —a, in the best one. 
Thus, from the microscopic point of view the two energy 
terms are cooperative. If hydrogen bonding is switched 
off, a = 0, we are back to the usual case of isotropic pair- 
wise attraction considered in Q , which yields a compact 
groundstate with no secondary structures. In the other 
limit of no hydrophobic interaction, a — oo, the ground- 
state has already been shown to be a long straight helix, 
when m > 6 10] . Since the ground state differs sig- 
nificantly in the two limiting cases, one should actually 
expect a non-trivial competition between the two energy 
terms for intermediate values of a, despite the micro- 
scopic cooperativity. This competition is induced at a 
global macroscopic level as a consequence of chain con- 
nectivity and excluded volume constraints [l2j . In this 
Letter, we will focus on its thermodynamic implications. 

Our results qualitatively agree with previous work on 
an analogous lattice model, where hydrogen bonding was 
mimicked via the introduction of rotating spins 131. We 
remark, nonetheless, that the extension of such results to 
our off-lattice model is higly non trivial, since the geo- 
metrical order implicit in the lattice structure could mask 
or enhance artificially secondary structure formation. As 
an example, while isotropic compaction of a homopoly- 
mer chain on a cubic lattice is sufficient to produce some 
amount of secondary structure Il4l , this is not true for 
an off-lattice homopolymer [Isl llfij . 

Our aim is to determine, by means of Monte Carlo sim- 
ulations, the density of states p (E) of a polymer chain 
with the Hamiltonian (jSJ), so that the partition function 
Zn (T) of a A-bead chain at reduced temperature T can 
be easily reconstructed: Zjy (T) = J2e P (E) ex P [~E/T\. 
We have employed a set of standard moves currently used 
in simulations of polymer chain; pivot, crankshaft, and 
reptation moves |l7j. In order to avoid trapping in lo- 
cal energy minima, we have employed a novel simulation 
method, based on generalized ensemble techniques (lif . 
The key notion, using generalized ensembles, is that a 
proper reweighting of temperature as a function of energy 
should allow the chain to escape from such energy min- 
ima 01 • The method lends itself in a natural way to be 
formulated within an iterative convergence scheme, and 
the possibility of properly employing the statistical infor- 
mation from more different steps of such scheme greatly 
increases its effectiveness 



Nevertheless, the presence of frustration provides an 
inherent limitation to such method, since it is based on 
the knowledge of local properties of the phase space, 
and the competition between different energy terms re- 
sults in different regions of the phase space sharing the 
same total energy but having different local densities of 
states. This turns out to be the case within our model, 
causing a very slow convergence to equilibrium. There- 
fore, we have introduced a finer, two-dimensional repre- 
sentation of the full multi-dimensional phase-space, by 
identifying a conformation through both its hydropho- 
bic energy Ehp and hydrogen bond energy Ehb- Our 
simulation method can be easily adapted in order to 
compute the density of states p (Eh p , Ehb) as a func- 
tion of both energy terms. Details on the employed 
reweighting scheme will be published elsewhere |20|. In 
this way, not only convergence to equilibrium is more 
easily obtained, but the partition function Zn (T, a) — 
Y,E hp ,E hb P( E h P , Ehb) exp [-{E hp + aE hb )/T], and hence 
any other relevant thermodynamic quantity, can now be 
reconstructed for any given value of the relative strength 
a between the two competing energy terms. The effec- 
tiveness of this sampling strategy shows that ther two 
energy terms serve as relevant order parameters. 

The main result of this work, obtained for a chain with 
A = 17 beads, is shown in Fig. 1. The logarithmic 
density of states S (E hp , E hb ) = In [p (E hp , E hb )], the mi- 
crocanonical entropy, is represented as a surface plot in 
the employed two-dimensional representation of the con- 
formational space. Effective free energy landscapes can 
easily be reconstructed within the same representation: 

Fa (T, Ehp, Ehb) = (Ehp + aEhb)/T — S (Eh p , Ehb) , (4) 

where the free energy F a (T, Eh p , Ehb) is given in dimen- 
sionless units. The surface plot in Fig. 1 can thus be 
interpreted as the opposite of the free energy landscape 
at infinite temperature, so that free energy valleys are 
seen as entropic ridges. In the phase space region with 
the lowest values of both energy terms, the entropy sur- 
face exhibits a rich yet regular structure, which is going 
to play a crucial role in determining the thermodynamic 
properties at low temperature. Three different ridges, 
separated by non-convex regions of the entropy surface 
corresponding to free energy barriers, can be identified. 

The properties of the conformation ensembles popu- 
lating such entropic ridges, or free energy valleys, can be 
readily identified by computing several order parameters, 
which measure the compactness degree and the amount 
of secondary structure content. Compactness is usually 
measured by means of the square gyration radius: 

N 

^ = £(rl-f cro ) 2 AV, (5) 

i=l 
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FIG. 1: Entropy surface plot S(Eh P ,Ehb) in the (Eh P ,Ehb) 
plane. Secondary structure content order parameters are 
shown in color scale, red for helices < E^ >, green for sheets 
< E(, s > + < E as >. The lighter the color the higher the 
order parameter. The yellow lines on the entropy surface 
show the average hydrophobic and hydrogen bond energies 
parametrized as a function of temperature for a = 0, 1,3, 4. 
The dashed portions of the curves refer to a first order tran- 
sition, which is identified by looking at the free energy contor 
plot, as in Fig. 3, and simply connect the two competing free 
energy minima. Typical conformations populating relevant 
entropic ridges are also shown. 

where f cm = J^iLi fii/N is the center of mass vector. As 
for secondary structures, we define the helical content of 
a conformation as: 

5 

X h = [(Vi-U-i + V i,i + V i+w) /3] m , (6) 

j— i=3 

and the parallel and antiparallel sheet content similarly: 

z PS = [(Yi-ij-i+Vij+Vi+ij+ii/sr > ( 7 ) 

j-i>6 

E QS = [(Yi-ij+i + Vu + Vw-ii/ST , (8) 

j— »>5 

where Vi t j — Vhb ifi — fj,Ui,Uj), (0 < Vij < 1) measures 
to what extent a hydrogen bond is formed between beads 
i and j. Each term in the above sums JBJ , (7J) , JSJl , again 
between and 1, measures to what extent the pair 
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FIG. 2: Specific heat per monomer C, black line, and mean 
square gyration radius < B? g >, blue line, as a function of 
reduced temperature T in logarithmic scale, in the a — 3 
case. In the inset, the specific heat per monomer, black line, 
is again shown together with the secondary structure order 
parameters < Eh >, red line, and < Eb s > + < E a3 >, green 
line. Dotted lines show the same quantities as computed from 
a second independent simulation. 



can be considered the center of a local portion of a given 
secondary structure. Within our definitions a simple hy- 
drophobic contact, which is not a good hydrogen bond, 
does not contribute to secondary structure counting. 

As is shown in Fig. 1, the ridges in the entropy surface 
are associated, with increasing number of hydrophobic 
contacts and decreasing number of hydrogen bonds, to 
helices with 4 beads per turn, to helices with 5 beads 
per turn, and to sheet-like conformations, respectively. 
The gyration radius decreases accordingly 20] , since long 
straight helices are extended objects. Whereas helices 
(a) and (b) are indeed representative of the two helical 
ridges, the sheet-like conformation (c) is just one among 
many possible different representatives. In this region 
we expect the occurrance of many different free energy 
minima, possibly giving rise to glassy behavior. Such 
frustration is of course not resolved within our bivariate 
parametrization. 

We now discuss in detail the case a = 3. The mean 
square gyration radius < R 2 g > , and the specific heat per 
monomer C = T~ 2 (< H 2 a > - < H a > 2 ) /N are shown 
in Fig. 2. The behavior of the secondary structure or- 
der parameters, < > and < E ps > + < S os >, is 
shown in the inset. The curves resulting from two differ- 
ent independent simulations are shown; the accuracy is 
quite good down to temperatures as low as 0.1. The spe- 
cific heat curve exhibits, with decreasing temperature, 
one shoulder and two higher and sharper peaks. Spe- 
cific heat peaks are usually related to a phase transition, 
but care should of course be taken in gene ralizing results 
from such a small system (see e.g. 21] for a detailed 
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discussion of problems arising in finite-size scaling of 6- 
collapse). For a finite-size analysis of such transitions, we 
show the free energy contour plots at the corresponding 
temperatures in Fig. 3. 

As signalled by the decrease of the gyration radius, 
the specific heat shoulder is related to the collapse of the 
chain from the swollen high temperature phase. The first 
peak close to the shoulder corresponds to a sharp increase 
of the gyration radius, and is related to the formation of 
many hydrogen bonds and to the appearance of helical 
structure, whereas the hydrophobic energy is almost not 
changed. The free energy contour plot clearly shows the 
existence of two competing minima, so this globule-to- 
helix transition is first order. The second peak is marked 
by a sharp decrease in the gyration radius, and is re- 
lated to the breaking of helices and to the formation of 
sheet-like structures, as (c) in Fig. 1. Sheet-like confor- 
mation are compact objects, having less hydrogen bonds 
but more hydrophobic contacts than helices. This last 
transition is again first order, as shown in the free en- 
ergy contour plot. Three different minima are present, 
originating from the entropic ridges identified in Fig. 1, 
but helices with 5 beads per turn (b) never get efficiently 
populated, since they suffer competition from either side. 

The very structure of the specific heat, one shoulder 
and then two peaks with decreasing temperature, is sim- 
ilar to what is found in the usual a = case, even if 
the thermodynamically stable phases are completely dif- 
ferent. It is tempting to interpret our results within the 
same overall framework proposed in [4J, that is to say 
with decreasing temperature the chain first undergoes a 
gas-to-liquid collapse, then a first order liquid-to-solid 
transition, and finally a solid-to-solid transition which 
is again first order (in the absence of hydrogen bonding, 
the last transition is a continuous polymorphic transition 
0). In our model, the possibility of hydrogen bonding 
simply acts in selecting helices and sheets among all pos- 
sible solid crystalline conformations. 

In Fig. 1 we have also summarized the different ther- 
modynamic static properties of the polymer chain when 
a is varied. The yellow lines can be thought of as dynam- 
ical trajectories only in the infinitely slow cooling case. 
Actual dynamics does not take place within the effective 
free energy landscape J3Jl , since kinetic barriers in the full 
phase space are smoothed over by the coarse-graining of 
our representation. This is most likely the case for the 
helix-to-sheet transition, where we expect the underlying 
energy landscape to be much more roughed. 

All trajectories in Fig. 1 start from a common point at 
infinite temperature, but then explore different regions of 
the phase space, according to different strengths of hydro- 
gen bonding. Nevertheless a collapse transition, related 
to the shoulder in the specific heat, is common to all a 
values, and moreover takes place in all cases for simi- 
lar values of the competing energies. At lower temper- 
atures, hydrogen bond strength greatly affects the ther- 
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FIG. 3: Contour plots at different temperatures, in the a = 3 
case, of the effective free energy F a (T, Eh P , Ehb), eq. ||1J, 
in the (Eh p ,Ehb) plane. The temperatures T = 0.46, T = 
0.14, correspond to the specific heat peaks seen in Fig. 2. 
The spacing between consecutive levels in each contour plot 
is unitary, and corresponds to a difference of ksT in physical 
units. The darker the color, the higher the free energy value. 
Letters refer to entropic ridges and conformations in Fig. 1. 



modynamic behavior. When a = 0, the conformational 
ensemble populated at low temperature does not include, 
as expected, the regions of the phase space with a high 
content of secondary structure. As in Q, two further 
transitions are present, liquid-to-solid and solid-to-solid. 
If a = 1, after the last solid-to-solid transition the sheet- 
like region of the phase space becomes efficiently sampled 
at low temperatures. If a = 3, as we have already seen, 
the chain first undergoes a transition from the liquid glob- 
ule phase to the helical region and then to the sheet-like 
region. Both transitions are first order. Finally if a = 4, 
hydrogen bonds are strong enough to produce a helical 
ground state, as in the a — oo case |lOj. 

Note that only the marginal border of the 'green' sheet- 
like region is thermodynamically relevant at low temper- 
atures (see also the small increase of the order parameter 
in Fig. 2). It is believed that a-helices are more likely 
to be formed by residues with small side chain groups, 
whereas the loss in conformational entropy suffered by 
bigger side chain groups, when arrenged in helical con- 
formation, favors the formation of /3-sheets j^. This 
general picture is consistent with our results. In fact, no 
side groups are present and helices are indeed entropi- 
cally favored, since they sit on the top of a ridge in the 
entropy surface, whereas sheet-like conformations do not. 

To summarize, we have introduced a simple model for 
an off-lattice self-avoiding polymer chain with two com- 
peting attractive inteactions, isotropic and directional- 
ized. By means of Monte Carlo simulations we have 
determined the density of states of the chain within a 
two-dimensional representation of the phase space, and 
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hence the phase diagram for different values of the rel- 
ative strength of the two competing energies. If the di- 
rcctionalized interaction is strong enough, different con- 
formational ensembles compete closely with each other 
at low temperature, which have peculiar proteinlike fea- 
tures, such as helices and sheets. 

It is a pleasure to acknowledge enlightening discussions 
with A. Maritan, C. Rischel, K. Sneppen and G. Tiana. 
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